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COMPUTER PROGRAMS FOR CALCULATION 
OF THERMODYNAMIC FUNCTIONS OF MIXING 
IN CRYSTALLINE SOLUTIONS 


I. INTRODUCTION 

Most of the important rock-forming minerals are crystalline solutions of 
two or more components. Therefore, it is necessary that mineralogists and 
petrologists become more familiar with the thermodynamic behavior of crystal- 
line solutions. Unfortunately the experimental data on the silicate solutions is 
meagre and quantitative calculations for many important minerals are not pos- 
sible at present. However, a semi -quantitative study of the data available from 
phase-diagrams and natural mineral assemblages may often be suitably used for 
a better understanding of the experimental and natural assemblages. The pro- 
grams described here are useful in various calculations for the thermodynamic 
functions of mixing and the activity-composition relations in minerals. These 
programs may be particularly useful to graduate students who may want to fa- 
miliarize themselves with thermodynamic behavior of solutions by computing 
various real or hypothetical problems. The thermodynamic equations used here 
are taken from Guggenheim (1952, 1967), Prigogine and Defary (1954), King 
(1969), and Saxena (1972). 

All of the programs used by Saxena (1972) are discussed below. The equa- 
tion numbers used by Saxena are given in brackets, [ ]. 
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PRECEDING PAGE BLANK NOT FILMED 


II. PROGRAM BETA 
A. Purpose 


This program may be used to solve the equation 

'i + <p 2A (p - i)l 


In \p 


l A 


zqi 

+ — — In 


z Ql 

= In « Pib + — In 


01A 0 3 + 1) 
1 + 02 B (P ~ 0 


01 B 03' + 1) 


[ 4 . 18 ] 


( 1 ) 


where \p 1A and t|/ 1B are the mole fractions of component 1 in A and B coexisting 
phases, and 0 2 are constant fractions defined by 


, x iqi , ^ 

1 x^! + x 2 q 2 2 Xjqj + x 2 q 2 ’ 

and (3 and p' are for A and B phases, respectively, and are given by the relation 

p = {1 + 40! 0 2 [exp (2W/ZRT) - 1]} * ; (3) 

q, and q 2 are constant factors and for very similar components, such as Fe 2 + 
and Mg 2+ may be taken as unity. Guggenheim (1944) considered zq , as the num- 
ber of sites which are neighbors of a molecule of type represented by component 
1. 


The notations \p lA and \J/ lB correspond to and x® and x 2A and \p 2B cor- 
respond to i //“ and x® in Saxena (1972). 


B. Numerical Method 


Setting y = 2W/ZRT, let 


zq x 

f(y) = — in 


~(P'(y) + 1) (P(y) + 1-2 <p 2A ) 
(P(y) + D(P'(y) + 1 - 20 2b ) 


+ In 




(4) 


Then the problem of finding y* such that (1) is satisfied becomes the problem of 
finding y* such that in (4) f(y*) = 0. 
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The method of solution is of bounding the zero, y*, above and below by y 
&.y 2 such that after the ith iteration 

ly i° - y? } l = (y ( i 0) - y 2 0) )/ 21 

where y* 0 * and y ( 2 0) are the initial bounds input to the program. Note: the 
assumption, 

y(0 < y* < y(*> 

is equivalent to 

fCy^ 0 ) f(y 2 1} ) < 0; y* (0) 

has initial value 


(y j 0) + y ( 2 0) )/2. 

At each iteration f(y*(')) is evaluated. If |f(y*(>))| < e ( in this program e = 
10 - 4 ), then the zero is considered found with y* = y* (i) ; else if 

f(yk 1} ) f(y* (l) ) < 0 


then 


y ( K 


i) = 


yk i} . 


v (i 1} 

J mod (k , 2) 


= y* (i) and y* (i + 1} = 


y ( i‘ 


i) 


y? 


i) 


for k = 1 or 2, provided i does not exceed a predetermined maximum, in which 
case the search for the zero is considered a failure. 


For each set (z, q t , q 2 , \p 1A , \p ) of data the program prints the following 
information: z, q lf q 2 , ^ 1A , ^ 1B , y*^ 0 , f(y* (i >), i, where y* = y* (i) if zero 
found else y*(0 is the final estimate when the search failed. 


There are two cases where failure can occur: 

(1) f(y ( ! 0) ) f(y ( 2 0) ) > 0; 

that is, y t and y 2 did not bound y*; 

(2) y(°) - y(°) 

was too large. 
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When either occurs additional information is printed as an aid: 

(1) y*< j >, f(y* (i) ), P(y* (i) ), p’ (y* (i) ), e****) 
for each value i assumed. 


C . Notation used in Program BETA 


X1A 

X 1A 

X2A 

X 2 A 

XI B 

X 1B 

X2B 

X 2 B 

Z 

Z 



Ql 

0i 

Q2 

^2 

PlA 

01 A 

P2A 

02 A 

P1B 

Yl 

PP 

P2B 

Y2 

02b 

y ( 2 0) 

Y(l) 

y 

Y(2) 

y? 

BETA 

P (y) 

BETAP 

& (y) 

F 

f(y) 



EY 

e* 



X 

current value of (S(y) 


XP 

current value of j3' (y) 


FY 

current value of f(y) 


ITER 

i 




Y3 = y* (0) 
Y(3) = y* (i) 


D. Input to and Output from Program BETA 


Card 1 


Card 2 


column 1-5 (right-adjusted) 

NX: number of pairs of (X 1A , X 1B ) to be evaluated for a 
given z, q 15 q 2 . 

columns 1-5, 6-10, . . . , 76-80 (fixed point format) 
(X1A, X1B) up to 8 pairs per card. Card 2 format is 
repeated until the NX pairs of (X1A, X1B) are entered 
8 to a card, except possibly the last. 


Last Card : columns 1-5, 6-10, . . . , 21-25 (fixed point format) 
Z, Yl, Y2, Ql, Q2 respectively 


This set of cards constitutes a case. Multiple cases are permitted, each 
case stacked one behind the other. 
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Figure 1 shows a sample set of input to program BETA while Figure 2 
shows a sample set of output. 


.990 .025 .975 .035 .965 .067 .933 .120 .880 
.95 1.05 

Figure 1. Sample Input to Program BETA 



z 

X 1 A 

X 1 0 

Y 


F ( Y) 

ft 

01 

02 

4. 

0 •£ J 5 

0 .99 5 

55673E 

V 1 

v .331 47E-«5 

1 5 

<• .95 • 

1 .0 59 

4« 

r .c iu 

C . 99*7 

^ • 49 726E 

0 1 

*.'• 4 2915 E-V4 

15 

C .9 56 

1 .0 53 

4. 

(■ .0 25 

6.975 

.:«.42 282F 

V* 1 

— f;* 44 823 E— i. 4 

14 

C .953 

1 .U 53 

4. 

0 .0 35 

C . 965 

l • 39 86 9E 

V*1 

-v* 972 75 E— C 4 

14 

9 .9 59 

1 .953 

4* 

U .06 7 

0. 933 

0.35272E 

0 1 

87738E-04 

14 

6 .9 59 

1 .053 

4. 

. 1 2* '. 

V . 856 

318 >2E 

31 

C ® 2 1 9 35E— J 4 

12 

< .9 57i 

1 .353 


.005 .995 .010 
4.0 0.0 6.0 


Figure 2. Sample Output from Program BETA 
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E. Listing of Program BETA 


C PROGRAM BETA 

C P . A .GOMEL L A 

C CODE 641,1 

C GODDARD SPACE FLIGHT CENTER 

C GREENBELT, MARYLAND 20771 

BET A( EY)=SQRT< 1 . + 4 ,*X1A*X2A*< EY-1 . ) > 

BETAP ( EY)=SQRT( 1,+4.*X1H*X2B*!EY-1. ) > 

F ! X * XP ) =ALOG( ( XP* 1 , ) * ( X* 1 ,-2.*P2A )/( ( X* 1 . ) * I XP+ 1 .-2 . *P2B ) ) ) - 2 2 
1 + ALOGI X 1 A/X 1 B I ♦ ALOG ( P 1 B /P 1 A ) *Z2 

COMMON STACK ( 6 t2 1 I , Y , F Y , X , XP ,E Y , I T , I T ER 
INTEGER*4 OUT/6/, IN/5/ 

REAL*4 A1{50>,B1<50),Y!3),EY(3),X!3),XP!3),FY<3> ,FYP(3> 

100 READ! IN, 1 , END= 1 500 ) NX, ( A 1 ( I) , B 1 (I) , I =1 , NX ) 

1 FORMAT ( 1 5 /(16F5.0I) 

READ! IN,5)Z,Y1»Y2,01,02 

5 FORMAT l 1 6F 5 ,0 I 
Y3=.5*(Y1+Y2) 

WRITE! OUT, 3 ) 

3 FORMAT! '1*,T4,'Z*,T13,'X1A' , T 22 , • X lb • ♦ T40, • Y « , T 5 3 , • F ( Y ) • ,T63, •#*, 

1 T73, 1 01 ' ,T83 * *02 * I 

300 Z2=.5*Z *01 

DO 1200 1=1, NX 
X 1 A= A 1 ( I ) 

X1B=B1! I ) 

X2 A= 1 • -X 1 A 
X2B=1 ,-XlB 

P1A=X1A*01/! X1A*Q1+X2A*Q2) 

P2A=1 ,-PlA 

P1B=X1B*Q1/! X1B*Q1*X2B*Q2 ) 

P2B=1 ,-PlB 

I T ER= 1 

IT =0 

Y ( 1 ) = Y 1 

Y! 2 ) =Y2 

Y(3)=Y3 

DO 500 J=1 , 2 

EY! J)=EXP! Y( JI/Z2) 

X! J ) =BETA! E Y( J ) I 
XP! J )=BETAP( EY! J) > 

FY< J)=F!X! J),XP(J)J 
I T=IT+1 
CALL SAVE ( J ) 

IF! ABStFY! J) ),LE, .0001 1 GO TO LOOO 
500 CONTINUE 

600 IF! ITER. LE. 20) GO TO 625 

WRITE! OUT ,4 ) Z,X1A,X1B, STACK 

4 F0RMAT(F4.0,5X,F6.3,5X,F6.3,T94, ' NON-CONVERGENC E * / T94 , 'BETA', 

1 T 105 , * BETA-P ' ,T114, 'EXP! 2Y/( Z*01 ) ) * ,T13'0, 'ITER' / 

2 (31X,2E13.5,30X,3E13.5,I5) ) 

GO TO 1000 

1 « X=A-INVERSE*B : LEAST SQUARES COE FF IC I ENTS ' /3020. 8 ) 

6 FORMAT! T 30* »K' ,T70, 'M* ,T90 * ' N * ,T50, * LN ( K ) • /T 20 , D20 . 8 , 20 X , 

1 2D20.8/T4, • J* ,T15, 'XAB* ,T35, • XAA • , T55, • Y • , T75, 

2 'Y-CALC* ,T95,‘R'/( 15,5020.8 ) ) 

625 EY(3)=EXP(Y< 3)/Z2) 

X! 3 ) =BETA( E Y! 3 ) ) 

XP(3)=BETAP!EY! 3) ) 

FY!3)=F(X!3),XP!3) ) 

I T= IT* 1 
CALL SAVE! J ) 

J = 3 

IF(ABS(FY!J)) .LE • .0001 ) GO TO 1000 
J2 = 2 

DO 700 Jl=l,2 
SIGN=FY( J )*FY( J1 ) 

IF! SIGN.GT.O. ) GO TO 700 
Y( J2)=Y< J) 

FY! J2)=FY! J) 

Y< Y< JI*Y( Jl) ) 

I TER= ITER* 1 
GO TO 600 
700 J2= 1 

1000 WRITE(0UT,2)Z,X1A,X1B,Y! J),FY<J), I TER, 01 ,02 

2 F0RMAT(F4.0,5X,F6.3,5X,F6.3,5X,2E13.5,5X,I5,2F10.3) 

1200 CONTINUE 

GO TO 100 
1500 RETURN 
ENO 

SUBROUTINE SAVEIJ) 

COMMON STACK! 5,21) , I STACK! 21 >,YSTACKt 3,5), IT, ITER 

00 100 1*1,5 

STACK! I , IT ) =YSTACK ( J , I ) 

100 CONTINUE 

ISTACK! IT » * ITER 

RETURN 

ENO 


00000100 

00000200 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000825 

00000810 

00000820 

00000815 

00000850 

00000875 

00000885 

00000950 

00001000 

00001100 

00001200 

00001300 

00001400 

00001410 

00001420 

00001430 

00001440 

00001500 

00001505 

00001525 

00001550 

00001575 

00001600 

00001700 

00001800 

00001900 

00002000 

00002050 

00002100 

00002200 

00002300 

00002350 

00002351 

00002352 

00002353 

00002354 

00002355 

00029250 

00032500 

00032700 

00032800 

00002360 

00002365 

00002370 

00002380 

00002383 

00002385 

00002390 

00002395 

00002400 

00002500 

00002600 

00002700 

00002710 

00002720 

00002730 

00002740 

00002750 

00002760 

00002770 

00003300 

00003400 

00003700 

00003800 

00003900 

00003905 

00003910 

00003915 

00003920 

00003925 

00003930 

00003935 

00003940 

0086 CARDS 
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PRECEDING PAGE , BLANK NOT FILMED 


HI. PROGRAM GEGIGM 

A. Purpose 

Program GEGIGM calculates the free energy of mixing, activity coef- 
ficients, ideal free energy of mixing and excess free energy of mixing in binary 
solutions. 

B. Numerical Method 

The program computes the various thermodynamic functions by using 
the following equations: 

RT lnf A = X| [A 0 + (3X A - X B ) + A 2 (X A - X B )(5X A - X B )] 

[1.48] 

RT lnf B = X\ [A 0 - (3X fi - X A ) + A 2 (X B - X A )(5X B - X A >] 

[1.49] 


G em - X* X B [Aq + Aj (X A - X B ) + a 2 (X a - x B ) 2 ] 

[1.53] 

G im = RT[X A lnX A + X B lnX B ] 
g m = g im + g em 

where f is the activity coefficient, G EM excess free energy of mixing, G IM ideal 
free energy of mixing and G M is the total free energy of mixing. 


A2 : A 2 

GAXA : f A X A 
GBXB : f R X B 
GM : G M 


D. Input to and Output from Program GEGIGM 

Card 1 : Columns 1-7, 8-14, 15-21, 22-28, 29-35 (fixed-point format) 
T, AO, Al, A2, R, respectively 
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Card 2 : Columns 1-5 

NX: number of observations of XA on which to perform the 
calculations. 1 < NX < 50. 

Card 3 : Columns 1-5, 6-10 66-70 (fixed-point format) 

XA: up to 14 per card. Card 3 format is repeated until the 
NX observations are entered, 14 to a card, except for pos- 
sibly the last card. 

These cards constitute a case. Multiple case are permitted, each case 
stacked one behind the other. 

For each case the following information is printed: 

1. T, R, AO, Al, A2 

2. for each observation: 

XA, XB, GA, GB, GE, GAXA, GBXB, GM 

3. Plots: GAXA vs. XA 

GBXB vs. XB 
GE vs. XA 
GM vs. XA. 

Figure 3 shows a sample set of input to GEGIGM while Figure 4 shows a 
sample set of output. 


1273 . 390.0 - 2177 . 0.00 1.987 

19 

0.03 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.53 0.60 0.65 0.70 
0.75 0.80 0.85 0.90 0.93 

Figure 3. Sample Input Data for Program GEGIGM 
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r=i27i. 

A0 = 

=oo . a i=* * - a ?= 

0 . 

X A 

X P 

CAMM A- A 

RAM s.»4--5 

35 

.0 50 

o.o=o 

? . cc 7 

1.007 

1 7 = . 74? 

. 1 00 

0 .0 00 

0 

• 

o 

o 

1.025 

2.7 5 , 044 

.150 

0 .550 

1.5=4 

1.0=5 

707,77? 

,?00 

0.5 00 

1 .705 

1.004 

751 .70? 

.750 

0.7=0 

1.710 

1 . 1 75 

770. 050 

.100 

0.700 

1 .002 

1.157 

750. ">5 0 

,150 

0.5 = 0 

1.007 

1 . ?^>5 

751 . 055 

.400 

0.50 c 

0.042 

1. 7 = 7 

71 = . 005 

.450 

0.5=0 

o .oo-j 

1 . 7?4 

?74. 155 

,5 00 

0 .500 

0 .5 5 1 

1.7=4 

7? ?. 500 

,550 

0.450 

0 .“7 1 

1 . 770 

155. 704 

500 

0.400 

0.57? 

1.757 

100.104 

550 

0 . 7 = 0 

0.55? 

1.74? 

= 7 . P05 

•700 

0.700 

0 .550 

1 . ?07 

4. 0.7? 

,75 0 

0.7=0 

0.015 

l.oig 

- 7 7 . ? 1 q 

,500 

0 .200 

0 .040 

1 . 1 22 

-f 5.40? 

,550 

0.1=0 

0.05? 

1 .005 
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Figure 4d 
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E 


Listing of Program GEGIGM 


c 

C PROGRAM GEGIGM 

C P.A.COMELLA 

C CODE 641.1 

C GODDARD SPACE FLIGHT CENTER 

C GREENBELT, MARYLAND 20771 

C 

REAL *4 XA( 50) ,XB( 50) ,GE( 50 ) ,GA< 50) ,GB< 50) ,GAXA< 5 0 ) , 

1 GBXBl 50 ) ,GM< 50 ) 

LOG I CAL* 1 GRID! 161,101 ) 

I NT EGER *4 NSC ALE ( 5 ) 75*0/ , NHL /8/ , NSBH/6/ , NVL /« / , NSBV/ 1 0/ 

1 ,NHL1/11/,NHL2/16/ 

100 READ(5»1,END=1000) T,A0,A1,A2,R 

1 FORMAT! 7F10. 3 ) 

RT = R*T 

200 READ (5,2) NX, ( X A ( I ) , I = 1 , NX ) 

2 FORMAT! 1 5 / < 14F5.3 ) ) 

DO 300 1=1, NX 

XB ( I ) = 1 XA ( I ) 

GAL = XB( I ) * XB ( I )*( A0+ A 1 * ( 3 . *X A ( I ) — XB ( 1 ) )+A2*< XA! I ) — X B ( I ) ) 

1 *( 5.*XA< I ) — XB ( I ) ) ) / R T 

GA ( I )=EXP( GAL ) 

GAXA! I ) =GA ( I ) *X A ( I ) 

GBL = X A ( I ) * X A ( I )*( A0-A1*! 3.*XB( I ) -X A ( I ) )+A2*( XB! I ) -X A ( I ) ) 

1 *( 5.*XB< I ) -X A ( I ) ) )/RT 

GE ( I ) =X A ( I ) *XB ( I )*( A0+! XA! I ) -XB ( 1 ) ) * ( A 1 + ( XA ( I ) -XB ( I ) ) * A2 ) ) 

GB( I ) = E X P ( GBL ) 

GBXB ( I ) =GB ( I ) *XB ( I ) 

G I = R T * ( XA! I )*ALOG(XA( I ) )+XB( I )*ALOG(XB( I) ) ) 

GM ( I ) =GI +GE ( I ) 

300 CONTINUE 

WR I TE ( 6, 3 ) T,A0,A1,A2,R 

3 FORMAT ( , 1T=',F5.0,5X, , A0='»F5.0,5X,'A1= , ,F5.0,5X,'A2= , ,F5.0, 

1 5X» 1 R= 1 , F6.4 ) 

WRITE! 6,4) ( XA! I ) ,XB( I ) ,GA( I ) ,GB( I ) ,GE( I ) ,GAXA( I ) ,GBXB( I ) ,GM( I ) , 
1 1=1, NX) 

4 FORMAT! ' XA 1 ,5X,' XB 1 , 5X , 1 GAMMA— A 1 * 5X , ' GAHHA-B 1 , 5X , 1 GE ', 

1 5X, 1 GAMMA! A ) *XA ', 5X, ' GAMMA < B )*XB ', 5X, • GM </ 

2 (F5.3,5X,F5.3,5X,F7.3,5X,F7.3,2X,F11.3,2X,F11.3,5X,F11.3, 

3 5X, F 1 1 . 3 ) ) 

CALL PL0T1 ( NSCALE ,NHL ,NSBH,NVL ,NSBV ) 

CALL PL0T2(GRID,1.,0.,1.,0.) 

CALL PL0T3! 'A 1 ,XA,GAXA,NX) 

WRITE (6, 5 ) 

5 FORMAT! ' 1 ' ) 

CALL PL0T3! 'B' ,XA,GBXB,NX) 

WR I TE ( 6 , 5 ) 

CALL PL0T4! 29, X A*GAMM A- A £ XB*GAMM A-B ' ) 

CALL PL OT 1! NSCALE, NHL l.NSBH, NVL, NSBV) 

CALL PL OT 2 ( GR I D , 1 . , 0 • , 450 • , — 1 00 . ) 

CALL PL0T3! • E ■ ,XA,GE,NX) 

WR I T E ( 6 , 5 ) 

CALL PL0T4! 11,' GE > ) 

CALL PL0T1(NSCALE,NHL2,NSBH, NVL, NSBV) 

CALL PL0T2! GRID, 1 . ,0. ,0. ,-1600. ) 

CALL PLOT 3 ( 1 M ' , X A , GM, NX ) 

WRITE! 6,5 ) 

CALL PL0T4! 11, • GM ' ) 

GO TO 100 
1000 RETURN 
END 


00000100 

00000200 

00000210 

00000220 

00000230 

00000300 

00000400 

00000500 

00000600 

00000700 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001900 

00002000 

00002100 

00002200 

00002300 

00002400 

00002500 

00002600 

00002700 

00002800 

00002900 

00003000 

00003100 

00003200 

00003300 

00003400 

00003500 

00003700 

00003800 

00003900 

00003925 

00003950 

00004000 

00004100 

00004200 

00004225 

00004250 

00004300 

00004400 

00004500 

00004600 

00004700 

00004800 
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IV. PROGRAM REGSOL1 
A. Purpose 


This program may be used to analyze the distribution of a component 
between two binary crystalline solutions which are now assumed to be "simple 
mixtures" (Guggenheim, 1967). 


B. Numerical Method 


The relation between X A , X^, co a , a/ , and k is given by 


X%(1 - X°) 
"(l - x{)(X a A ) 


W a W0 

= InK + (1- 2X“) (1 - 2X0) 

RT RT 

[3.8] 


(5) 


where X's are mole fractions of A and B in a and j3, W’s are "interchange" en- 
ergies and k, the equilibrium constant. 


Given (5) and a set of NX observations (X“ if X^ j), i = 1, 2 NX, the 

problem is to find the best estimates for 


w a 
k ’ RT 


W 0 

, and—, 


according to the method of least squares. Let 


W a 


M = — N = W^/RT, X t = 1, 


X, = 1 - 2 X® X 3 = 1 - 2 X{ 


Y = In 


a - xj) x0 


X A O - X i> 
k’ = In (k) 


(5) can be rewritten as 


y = k'X 1 + NX 2 - MX 3 _ (6) 

The set of original observations (X® if X^.) are now transformed into the se- 
quence of observations (X^ , X 2i , X 3i , Yj), i = 1, 2, . . . , NX which can be 
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used in (6) to obtain the coefficients, k\ N, M, in the same way as out- 
lined in PROGRAM MATRIX. 


C. Notation used in Program REGSOL1 


XAA 

XAB 

XAA2 

XAB2 

KO 

MO 

NO 

CAY 

Y 

YEST 

R 

KCALC 

RK 


xi 

X 2 

X 3 

k' 

M 

N 

k 

y 

y as calculated using the least squares coefficients, k', M, 
N 

y estimate - y 

k - calculated from (6) holding M, N constant 
R calculated - k 


CHISQ 


2 RVY 

l 


KCHI : 2 RK 2 /k 


D. Input to and Output from the Program REGSOL1 


Card 1 : column 1-5 (right-adjusted) NX: numbers of pairs (X^, Xg) 

Card 2 : columns 1-10, 11-20, . . . , 71-80 (X® , X®) up to 4 pairs 
per card. Card 2 format is repeated until the NX pairs 
(X®, X®) are entered 4 to a card, except possibly the last. 

These cards constitute a case. Multiple cases are permitted, each case 
stacked one behind the other. 

For each case the following information is printed; 

(1) The least squares matrix, A, by column, 

(2) The B vector (the solution X = A -1 B). 

(3) A -1 , X which contains the least squares coefficients. 
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(4) CAY, KO, MO, NO, (J, XAB(J), XAA(J), y(J), YEST(J), R(J), J - 1, NX). 

(5) CfflSQ. 

(6) (J, CAY, KCALC(J), RK(J) , J = 1, NX). 

• (7) KCHI. 

Figure 5 shows a sample set of input to REGSOL1 while Figure 6 shows a 
sample set of output. 


9 

0 . 0 ? 1 
0.258 
0.029 


0 . 341 
0.902 
0 . 539 


0.070 

0.341 


0.692 

0.899 


0.094 

0.533 


0.815 

0.907 


0.136 

0.033 


0.864 

0.475 


Figure 5. Sample Input to Program REGSOL1 



■MATRIX (BY COLUMN): 

LEAST SQUARES MATRIX 



0.900000000 

01 

0. 597 000 OOD 

0 1 

-O.3068COOOO 

0 1 

- 0 • 59 70 0 0 000 

01 

-0.49461 480 J 

01 

0 . 1 6830920D 

01 

- 0. 38680 0 0 00 

01 

-o. 168309200 

-0 1 

0 . 31 2986400 

01 

VECTOR 






- 0.292371 030 

02 

—0.205707 02 D 

02 

0.1211824 60 

02 

■INVERSE (BY COLUMN) 





0.210415660 

01 

0 .202540020 

01 

0.1511 1833D 

01 

-0.202540620 

01 

-0.21972 1 69Q 

01 

-0. 1 32161 090 

01 

0.151 11 8330 

01 

0. 132lol 090 

01 

0. 1 47637860 

01 

- I NVER S~.* 0 1 LEAST SQUARES COErFICIt 

MTS 



- 0. 15408793'"* 

01 

0 • 199 45 51 3D 

01 

0.394962610 

00 


to 

o 


Y=LN( XAQ*( 1-X4A)/(XAA*( l-XAd)> 

Y-CALC=LN( K >-M* ( 1-2 «X AES >+N*< 1 -2* AAA) 

R = YEST-Y 

K LN( < ) M N 



0. 

214192680 00 

-0 • 

154087930 01 

0 . 

1^>9455180 01 

0 

• 89496261 D 00 


J 

XAR 

X A A 


Y 


Y-CALC 


R 


1 

0.21 0000000-01 

0. 3*1000000 

00 

-0 . 31 831 681 0 

01 

-0.316706180 

01 

0. 1 61 063880-01 


2 

0. 700000000-0 1 

0. 69200000D 

00 

-0. 33961 755? 

01 

-0. 359985940 

01 

-0.203683920 00 


.3 

0. 940000000-01 

0. 815000 0 00 

00 

-0. 374857680 

01 

-0 . 3724281 6D 

01 

0.242950520-01 


4 

0. 1 3600 0000 00 

0. 364000000 

00 

-0 • 369733580 

01 

-0 . 364444570 

01 

0. 5339001 8D-01 


5 

0.258000009 00 

0. 902 000 00D 

00 

-0.327603670 

01 

-0. 322579230 

01 

0.502444230-01 


6 

0.34100000? 00 

0. 699000000 

00 

-0 . 284500360 

0 1 

-0 .283932690 

01 

-0. 443233330-01 


7 

0.5 3300000? 00 

u. 90 7 00 0 0 CD 

00 

-0.214535000 

0 1 

-0 .213773840 

01 

0.761235850-02 


8 

0. 3300QC00D-01 

0. 47500000? 

00 

-0 . 327760750 

01 

-0. 33590425D 

01 

-0.814350220-01 


9 

0. 290000000-01 

0.5390 00 0 OC 

00 

-0 . 366734820 

01 

-0.3439554 ID 

01 

0.177794040 00 


CH I SQ= 

SUMM AT I UN ( R**2/Y 

)= -0.253566160-01 






K-CALC 

= <£XP( M 0 * X A 3 2 ) /EXP (N0*XA A2 ) ) *YEXP 








RK=K-CALC - KO 









J 




KO 


K-CALC 



RK 

1 




0.214192680 

00 

0.210770450 

00 


-0. 342223660-02 

2 




0 .214192680 

00 

0.262581 080 

00 


0.483884000-01 

3 




0 .2 14 19 268.0 

00 

0.209051570 

00 


-0.5141 l 1760-02 

4 




0.214192680 

00 

0 • 20 305685D 

00 


-0. I 11 35835D-0 l 

5 




0.214192680 

00 

0.203696590 

00 


-0.104960950-0 1 

6 




0 . 2 14192680 

00 

0.223899960 

00 


0.970727380-0 2 

7 




0.214192680 

00 

0.212568360 

00 


-0 .16243212D-02 

8 




0.21 4 192 680 

00 

0.232365370 

00 


0. 1817269 10-0 1 

9 




0.214192680 

00 

0 • 179303870 

00 


— 0 .348888 14D— 0 1 


CHI SQ = S UM M A T I DN( RK* “2/KO )= 0.1 90797660-0 1 


Figure 6. Sample Output from Program REGS0L2 



o o o o 


E. Listing of Program REGSOL1 


c 

200 

2 

300 

3 

600 

700 

7 

1 


PROGRAM :REGS0L1 
P.A.COMELLA 
CODE 641.1 

GODDARD SPACE PLIGHT CENTER 
IMPLICIT REAL*8 <A-H,0-Z) 

REAL *8 K(10),M<10>,N(10),XAA(100),XAB(100),XAA2<100), 

1 XAAQ! 100) ,XAB2< 100) ,XABQ< 100 ) , XABA t 100) , Y ( 100 ) , 

2 YEST(IOO) ,KO,MO,NO,R< 100 ) , A< 3,3 ) ,B< 3 ) ,YEXP(100) 

4 tKCALCI 100 ) ,RK( 100 ) ,KCHI 

COMMON K0,M0,N0,XAB,XAA,XAB2,XAA2,NX 
I OU T = 6 
I N= 5 

INPUT 

READ! IN, 2,ENI)=1600) NX,(XAB(I),XAA(I),I=1,NX) 

FORMAT ( 15/1 HF10.3 ) ) 

DO 300 1=1, NX 
X A A 2 ( I ) =1 .DO-2 .DO*XAA< I ) 

X A AO ( I ) = XAA2 ( I )*XAA2< I ) 

XA6 2 ( I >=1.00-2.D0*XAB< I ) 

X ABO ( I ) = X AB2 ( I ) 4= X A B 2 ( I ) 

XABA ( I ) =XAB2 ( I )*XAA2< I ) 

Y EXP ( I ) = < (XAB< I )#< 1 . DO-X AA ( I ) ) ) / ( XA A ( I ) * ( 1.D0-XAB1 I ) ) ) ) 
Y( I )=DLOG( YEXPI I ) ) 

CONTINUE 
JGO= 1 

WRITE! I OUT, 3) 

FORMAT! '1 PROGRAM REGSOL 1 ' ) 

DO 600 IM=1,3 

B ( I M ) =0 .DO 

DO 600 JM= I M , 3 

A! I M , JM ) =0 .DO 

00 700 1=1, NX 

A! 1 ,2 ) = A( 1 , 2 ) + XAB2 ( I ) 

A! 1 , 3 ) = A ( 1 , 3 ) +X AA2! I ) 

A(2,2)=A(2,2)+XAB0( I ) 

A ( 2,3 ) =A( 2, 3 ) +XABA ( I ) 

A ( 3 , 3 ) = A ( 3,3) +X AAO ( I ) 

B( 1 )=B( 1 )+Y( I ) 

B< 2 ) = H( 2 )+Y< I )*XAB2( I ) 

B ( 3 ) = B ( 3 I+XAA2! I ) * Y ( I ) 

CONTINUE 
A! 1,1 ) =N X 
A! 2 ,1 )=A( 1 ,2 ) 

A! 1 , 2 ) =- A ( 1,2) 

A! 2 ,2 )=-A< 2,2) 

A ( 3 , 1 ) = A ( 1,3) 

A ( 3 ,2 ) =-A( 2,3) 

WRITE! I OU T , 7 ) A , B 

FORMAT! 'OA-MATRIX (BY COLUMN ): LEAST SQUARES MATRIX' /3D20. 
1 3D20.8/' B-VECTOR ' /3U20.8 ) 

CALL MAT IMV! A ,3 ,B , 1 , DETERM ) 

WRITE! I OUT , 1 ) A , B 

FORMAT! 'OA-INVERSE (BY COLUMN ) 1 /3D20 . 8/ 3D20 . 8/ 3020 . 8 / 

K 0 = B ( 1 ) 

M0 = B ( 2 ) 

NO=B ( 3) 

CAY = DEXP(K0 ) 

CHI SO=O.DO 
KCHI=O.DO 
DO 800 J = 1 , N X 


00022600 

00022700 

00022800 

00022900 

00023010 

00023100 

00023300 

00023400 

00023800 
00023900 
00024000 
00024100 
00024200 
00024300 
00024400 
00024500 
00024600 
00024700 
00024800 
00025700 
00025800 
00025900 
00026500 
00026600 
00026700 
00026800 
00026900 
00027000 
00027100 
00027200 
00027300 
00027400 
00027500 
00027600 
00027700 
00027800 
00027900 
00028000 
00028100 
00028200 
00028250 
00028300 
00028800 
8/3D20. 8/00028900 
00029000 
00029100 
00029200 
00029225 
00029300 
00029400 
00029500 
00029600 
00029800 
00029850 
00029900 
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YEST(J) = K0-M0*XAB2( J )+N0*XAA2< J) 

R(J)=YEST(J) -Y(J) 

CHISO= R ( J ) **2 / YIJI+CH1S0 

KCALC! J) =DEXP( MO*XAB2( J >-NO#XAA21 J > )*YEXP( J ) 

RK( J ) =KCALC ( J ) -CAY 
KCHI=KCHI+RK< J)**2/CAY 
800 CONTINUE 

WRITE! I OUT, 10) 

10 FORMAT!////' Y=LN(XAB*< 1 -X A A ) / ( X AA* ( 1-XAB)) •/• Y-CALC=LN( K ) -M* • 

1 T18,'!1-2*XAB)+N*(1-2*XAA)'/' R=YEST-Y*> 

850 WRITE! I0UT,6) C A Y , KO , MO , NO , ( J , X AB ( J ) , 

1 X A A ( J ) , Y< J) , YEST! J ) ,R( J ) ,J=1,NX) 

WRITE! I OUT , 9 ) CHISO 

9 FORMAT! ' CHISO=SUMMATION< R**2/Y )=',D20.8) 

WRITE! I OUT, 1 1 ) ( J ,CAY ,KCALC( J ) , RK ( J ) , J = 1 , NX ) 

11 FORMAT!/' K-CALC=(EXP(M0*XAB2)/EXP!N0*XAA2) >*YEXP'/ 

1 ' RK=K— CALC - KO'/ 

2 T4,'J',T55,'K0',T75, ' K-CALC ' , T1 1 5, 'RK ' / 

3 < I5,40X,2D20.8,20X,D20.8) ) 

WRITE! I OUT , 12 ) KCH I 

12 FORMAT!' CHISQ=SUMMATI0N!RK**2/K0)=',D20.8> 

1000 CONTINUE 

1200 CONTINUE 
1500 CONTINUE 

GO TO 200 
1600 RETURN 
ENO 

SUBROUTINE MATINV! A, N,B,M, DETERM ) 

C MATINV IS A VERSION OF THE SHARE SUBROUTINE OF THE SAME NAME 

IMPLICIT R E A L * 8 (A-H,0-Z) 

RE AL*8 A(N,N) , BIN, M) , PIVOT! 10) 

INTEGER*^ IPIVOT! 10 ), INDEX! 10,2 ) 

EQUIVALENCE ( IROW,JROW) ,( I COLUM, JCOLUM ) , ( AM AX , T , SWAP ) 

DETERM= 1 . DO 
DO 20 J = 1 , N 
20 IPIVOT (J)=0 
DO 550 1=1, N 
AMAX=O.DO 
DO 105 J= 1 , N 

IF! IPIVOT! J) .EQ. 1 ) GO TO 105 
DO 100 K = 1 , N 

IF! IPIVOT(K)-l) 80,100,740 
80 IF! DABS! AMAX ) .GE.DABS! A! J ,K ) ) ) GO TO 100 
I ROW = J 
I C OLUM =K 
AMAX = A ( J , K ) 

100 CONTINUE 
105 CONTINUE 

IPIVOT! ICOLUM)=IPIVOT( ICOLUM ) + l 
IF! IROW.EQ. ICOLUM) GO TO 260 
DETERM=-DETERM 
DO 200 L = 1 , N 
SWAP = A( I ROW , L ) 

A! I ROW , L ) = A ( ICOLUM, L) 

200 A! ICOLUM, L ) =SWAP 

IF(M.LE.O) GO TO 260 
DO 250 L=1,M 
SWAP = B ( I ROW , L ) 

B ( I ROW , L ) = SW AP 
250 B( ICOLUM, L )=SWAP 


00030000 

00030100 

00030400 

00030410 

00030420 

00030430 

00030500 

00030600 

00030700 

00030800 

00031000 

00031100 

00031900 

00032000 

00032010 

00032020 

00032030 

00032040 

000320501 

00032060] 

00032070! 

00032300! 

00032400] 

00033000j 

00033100 

00033200 

00033300 

00009700 

00009800 

00009900 

00010000 

00010100 

00010200 

00010300 

00010400 

00010500 

00010600 

00010700 

00010800 

00010900 

00011000 

00011100 

00011200 

00011300 

00011400 

00011500 

00011600 

00011700 

00011800 

00011900 

00012000 

00012100 

00012200 

00012300 

00012400 

00012500 

00012600 

00012700 

00012800 
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260 INDEX! I ,1 )=IR0W 

INDEX! 1 , 2 ) = I COLUM 
PIVOT! I ) = A ( I COLUM, I COLUM) 

DETERM = DE TERM* PIVOT! I ) 

A! ICOLUM, IC0LUM)=1.D0 
DO 350 L = 1 , N 

350 A! ICOLUM, L ) =A( ICOLUM, L ) /PI VOT ( I ) 
IF(M.LE.O) GO TO 380 
DO 370 L = 1 , M 

B( ICOLUM, L)=B< ICOLUM, L) /PI VOT! I! 
370 CONTINUE 
380 DO 550 L1=1,N 

IF! LI. EQ. ICOLUM) GO TO 550 
T= A ( LI, ICOLUM) 

A! LI, ICOLUM ) =0.00 
DO A 50 L = 1,N 

A50 A! L 1 ,L ) =A( L 1 ,L )-A( ICOLUM, L)*T 
IF(M.LE.O) GO TO 550 
DO 500 L = 1 , M 

500 B(L1,L)=B(L1,L) — B( ICOIUM. L)*T 

550 CONTINUE 

DO 710 1=1, N 
L = N+ 1 - I 

I F( INDEX! L, 1 ) .EQ. INDEX! L,2 ) ) GO TO 
JROW= INDEX! L , 1 ) 

JCOLUM= INDEX! L , 2 ) 

DO 705 K = 1 , N 
SWAP=A(K, JROW) 
A(K,JROW)=A(K,JCOLUM) 
A(K,JCOLUM)=SWAP 
705 CONTINUE 
710 CONTINUE 
7A0 RETURN 
END 


00012900 
00013000 
00013100 
00013200 
00013300 
00013A00 
00013500 
00013600 
00013700 
00013800 
00013900 
0001A000 
000 1 A 1 00 
0001A200 
000 1 A 300 
0001AA00 
000 1 A 500 
0001A600 
000 1 A 7 00 
0001A800 
0001A900 
00015000 
00015100 

710 00015200 

00015300 

00015A00 

00015500 

00015600 

00015700 

00015800 

00,015900 

00016000 

00016100 

00016200 
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PRECEDING page blank NOT filmed 


V. PROGRAM REGSOL2 
A. Purpose 

If the data on K, W a ,/RT, W0 /RT are available, we may calculate X“ 
or X^, (given one or the other) in (5) and plot these on a Roozeboom figure. 

This provides us with a distribution curve or isotherm representing the distribu- 
tion of a component between two binary solutions 


B. Numerical Method 

REGSOL2 assumes that in (5) K, W a /RT, W0/RT and X^ are given, the 
problem then being to find X“ . To accomplish this (5) is transformed as 
follows : 


g(X“) = (1 - X£)exp(-Nn - 2X“» - f(Xj)X£ (7) 


where 


M = W a /RT, N = W0/RT, 


f(X^) = 


exp 


In f k 


(i - 

A 


x^ 


-M(l -2X5) 


The Newton-Raphson method is then applied to (7) with 

x? 0) = x{/(x{ + k (1 - x{)) 


the zeroth estimate of X^ . 

Each subsequent estimate is given by 


where 


X A° + 1} = X*(0 


g(Xf)) 
g' (X®* 1 )) 


( 8 ) 


g' (X£> = (2N (1 - X£) - 1) exp (— N ( 1 - 2X)) -f(X^) 

Whenever |g(X®('))| < e, (e is set in the program at 10 -4 ), the zero is said to 
have been found. This method fails in that region of the X“ vs. X0 curve where 
the slope is parallel to the X“ axis. 
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C . Notation used in Program REGSOL2 

XAB 


XAAC 

C 

x 5 

f<x£) 

K 

K0 

M 

M0 

N 

NO 

ALK 

k - calc 6 using XAAC, XAB, M, N in (6) 

R 

K0 - ALK 

D. Input to and Output from Program REGSOL2 


Card 1 : column 1-10, 11-20, 21-30 (fixed point format) K0, MO, NO, 

respectively. 

Card 2 : column 1-5 (right-adjusted) NX: number of observations 

for which X^ is to be found. 

Card 3 : columns 1-10, 11-20, . . . , 71-80 (fixed point format) X^ 
up to 8 observations per card. Card 2 format is repeated 
until the NX values of x(* are entered 8 to a card, except 
possibly the last card. 

This set of cards constitutes a case. Multiple cases are permitted, each 
case stacked one behind the other. 

For each observation the following output is provided if the zero has been 
found: 

X^, k Q , k-calculated, R. 

For the entire case a graph of X® versus is given. 

When the Newton-Raphson method fails to find the zero the following informa- 
tion is given: 

X A< X A (0 > f < X A (i) ) 

for each iteration i. 

Figure 7 shows a sample set of input to program REGSOL2 while Figure 8 
shows a sample set of output. 
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0.2823 2.0000 1.3390 

16 

0.03 0.05 0.06 0.08 0.10 0.15 6.20 0.25 

0.30 0.40 0.50 0.60 0.70 0.80 0.90 0.95 

Figure 7. Sample Input to Program REGSOL2 


FOB XAB= 0.25000E OO ZERO WAS NOT FOUND. TRACE* OF ITERATIONS FOLLOWS: 
X Ft X » 

0.541 A4782E 00 0.34369034E 00 

0.65933743E 01 -0.68313744E 08 
0.6 24 33 30 OE 01 -0.250e0240E 08 
0.5 89 47 4 2 0 E 01 -0.92051840E 07 
0 . 554 7796 2 E 01 -0.33774710E 07 
0.52027168E 01 -0.12387390E C7 
0 . 4859773 6 E 01 -0.45410481E 06 
0.451 92976E 01 -0.16636613E 06 
0.4181 70 1 7E 01 — 0. 6090 264 8E 05 
0.38475027E 01 -0.222724S7E 05 
0.35173635E 01 -0.81342930E 04 
0.319214S3E 01 -0.29654358E 04 
0 . 2 872984 9E 01 -0.107838C1E 04 
0.25614128E 01 -0.39076123E 03 
0.22595224E 01 -0. 14085666E 03 
0 . 19702377E 01 -0.50367310E 02 
0. 16977329E 01 -0.1777E421E 02 
0.14481325E 01 -0.61282244E 01 
0.1 2306871E 01 -0.201S8e63E 01 
0 . 10592642E 01 -0.59501708E 00 

Figure 8a. Sample Output from Program REGSOL2 Trace When Zero Not Found 



K= 0.2823E 00 


0*20 00 E 01 

N = 

0. 1339E 01 





NF WTCN- 

paphson method: 

ITERAT ION 

0 3 







J 

X AB 


XAA 

XAA (CALC) 


K (CALC ) 


K( INPUT) 


K-K(CALC) 

1 

0.3000E-01 


0.2908E 

00 

O.2023E 

00 

0.2823E 

00 

0. 1 3UE-05 

2 

0. 5000E-01 


0.5Bft2= 

00 

0.2823F 

00 

0.2823E 

00 

-0.2980E-06 

3 

0 « 6000 E- 

•0 1 


0. APOOE 

00 

0.2822E 

00 

0.2823E 

00 

0.5084E-04 

4 

0.8000E-01 


0 . 7 75 7E 

00 

0.2823E 

00 

0.2823E 

00 

0. 184BE-05 

5 

0. 1000E 

00 


0 • A220E 

00 

0.2823E 

00 

0.2823F 

00 

0.4292F-0S 

6 

0. 1 500F 

00 


xxxx 

ZERO 

NOT FOUND 





7 

0.2000E 

00 


X XX X 

ZERO 

NOT POUND 





ft 

0.2 500E 

00 


xxxx 

7FR0 

NOT FOUND 





9 

0 .3000E 

00 


0.0 102= 

00 

0.2822E 

00 

0.2823E 

00 

0.B500E-04 

10 

0 • 400 OE 

00 


0.9 1455 

00 

0.2823F 

00 

0.2823F 

00 

0 • 4548 F— 04 

11 

0 .5000E 

00 


0.9150E 

00 

0. 282 3F 

00 

0.2823F 

00 

0. 1B48E-05 

12 

0.6000E 

00 


0.9155= 

00 

0.2823E 

no 

0.28P3E 

00 

0.6139F-05 

13 

0.7000E 

00 


C.9195E 

00 

9.2823E 

00 

0.2823E 

00 

0.4089E-04 

14 

0. 800 OE 

00 


0.931 3E 

00 

0.2822E 

00 

0.2823F 

00 

0.7963F-04 

IS 

0 • 9000E 

00 


0.9562“ 

00 

0. 2823F 

00 

0.28235 

00 

0.3219F-05 

16 

0.9500F 

00 


0.9 755 = 

00 

0. 2823E 

00 

0.2823F 

00 

0. 1 848F-05 


Figure 8b. Sample Output from Program REGS0L2 
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1.000 


0. 500 


0.375 


0 . 250 ♦ 


0. 125 ♦ 


0.000 +- 


X 

X 

I 

I 

I 

I 

I 


I 

* X 

I 

X 

I 

1 

I 

I 

I 


I * 

X 

f 

I 

X 

I 

X 

I 

I 

* 

X 

I 

I 

I 

I * 

1 * 

* 

* I 

* I 


X 

X 

X 

I 

I 

X 

X 

1 

I 


I 

X 

I 

X 

I 

I 

t 

I 

X 


I 

I 

I 

I 

I 

I 

X 

I 

I 


X 

X 

I 

* X 

I 

I 

I 

I 

I 


X 

I 

I 

I 

I 

I 

X 

I 

I 


T 

I 

I * 

X 

I 

I 

I 

X 

I 


X 

1 

X 

I 

I 

I 

I 

I 

I 


X 

I 

I 

I 

X 

I 

I 

X 

I 


X 

X 

I * 

I 

I 

I 

I 

I 

I 


I 

X 

I 

I 

I 

I 

I 

X 

X 


X 

y 

I 

I 

I 

I 

I 

X 

I 


I 

X 

I 

I 

I 

I 

I 

I 

I 


I 

X 

I * 

I 

I 

1 

1 

X 

I 


X 

X 

z 

I 

I 

I 

t 

I 

T 


I 

I 

I 

I 

X 

X 

I 

I 

X 


X 

X 

X 

I 

I 

I 

X 

X 

I 


I 

X 

I 

1 

I 

I 

I 


X 


I 

X 

I 

X 

I 

X 

I 

I 

I 


I 

I 

I 

I 

I 

I 

I 

I 

I 


I 

I 

I 

I 

I 

X 

I 

I 

T 


* 

X 

X 

I 

X 

I 

I 

I 

I 


X 

I 

I 

X 

I 

I 

I 

I 

I 


I 

X 

X 

I 

I 

I 

I 

I 

I 


X 

I 

I 

I 

I 

I 

I 

I 

I 


X 

X 

X * 

I 

I 

I 

I 

I 

I 


I 

? 

I 

X 

I 

I 

X 

X 

X 


T 

I 

I 

I 

I 

I 

I 

I 

X 


X 

X 


I 

I 

X 

X 

X 

I 


I 

X 

X 

I 

I 

I 

I 

I 

I 


I 

X 

I 

I 

X 

X 

I 

I 

X 


I 

I 

I 

X 


I 

I 

X 

I 


I 

X 

I 

I 

I 

t 

I 

I 

I 


I 

I 

I 

I 

1 

1 

X 

X 

X 


I 

X 

I 

I 

I 

I 

I 

I 

I 


I 

I 

I 

( 

I 

I 

I 

I 

X 


I 

X 

I 

I 

I 

I 

I 

X 

X 



X 


K 0 = 


0.0 


0 . 282 30 E 00 


125 

XAB 


0.375 


0.750 


0.875 


1.000 


M0= 


0.20000E 01 


N0= 0.I3390E 01 


Figure 8c. Sample Output from Program REGS0L2 Plot of XAA vs XAB 
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E. Listing of Program REGSQL2 


C 

c 

c 

c 

c 

c 


c 

100 

2 

3 


1 

1000 


100 


4 


200 


500 


1 


7 

600 

8 

700 

3 


2 


5 


PROGRAM REGS0L2 
P.A.COMELLA 
CODE 641,1 

GODDARD SPACE FL 1 GHT ' CENT ER 
GREEN BELT , MARYLAND 20771 

REAL *4 K,M,N 

COMMON, K,M,N,X( 100) ,NX 
INPUT 

READ(5,2,END=1000) K,M,N 
FORMAT ( 3F10.3) 

READ! 5*3) NX , ( X! I ) , I =1 ,NX ) 

FORMAT ( I 5/ ! 8F 10. 3 ) ) 

CALL ESTMTE 
WR I TE ( 6* 1 ) 

FORMAT ( * 1 PROGRAM REGS0L2 * ) 

GO TO 100 
RETURN 
END 

SUBROUTINE ESTMTE 
IMPLICIT REAL*4 !A-H,K-Z> 

INTEGERS ITER /20/ , IOUT/6/ *FLAG 
1 ,NSCALE(5)/ 5*0/ * NHL/8/ » NSBH/6/ »NVL/8/,NSBV/10/ 

REAL*4 XAAC! 100) ,R( 100) »T0L/1. 0-4/ ,STACK(2,20) 

COMMON K0»M0,N0 ,XAB(100) »JX 
REAL*4 XAB21 100) ,RESID/0./ , X ALF / ' XXXX » / 

DIMENSION XXAB! 100) , ALK! 100) 

LOG 1CAL*1 GRID! 5200 ) 

FUN ( K, M, X, X2 ) = EXP(ALOG(K*( ( 1.00-X)/X) )-M*X2) 

GFUN ( N , X ) = ( 1.00-X)* EXP(-N*( 1.00-2 -00*X) )-C*X 

GPF ( N, X ) = ( 2. 00*N* ( 1.00-X 1-1.00)* EXP(-N*( 1 .00-2.00*X ) ) -C 

DO 500 J=1,JX 

XAB2( J )= 1.00-2. 00*XAB( J ) 

C=FUN(K0»M0»XAB( J) *XAB2( J) ) 

X = X A8 ( J)/< XAB( J)+KO*( l.-XABI J ) ) ) 

DO 100 JTER= 1 » ITER 
GX = GFUN( NOt X ) 

STACK! ltJTER)=X 

STACK! 2 * JTER ) =GX 

IF! ABS! GX) .LT.TOL) GO TO 200 

X=X-GX/GPF(NO,X) 

CONTINUE 
XAAC! J)-XALF 

WRITE! 1 OUT * 4 ) XAB ( J ) » ST ACK 

FORMAT! f OFOR XAB= ' t E 1 3. 5 t ' ZERO WAS NOT FOUND. TRACE' 

1 ' OF ITERATIONS FOLLOWS: ' /T12» 'X • ♦T28,‘FI X )' /!2E16. 8) ) 

GO TO 500 
XAAC! J )=X 

ALK ( J ) =ALOG( ( l.-X) *XAB ( J ) / ! ( 1 .-XAB ( J ) > *X > ) +MO*X AB2 ( J ) 

1 -NO*! 1.-2. *X) 

ALK! J)=EXP( ALK! J) ) 

R ( J ) =KO-ALK ( J ) 

CONTINUE 
RES I D =0 . 

WRITE! I OUT, 1 )K0, MO, N0» JTER 

FORMAT ( ' 1 PROGRAM REGSOL-2’/ 

1 T10, 'K=* ,E12.4,T30,'M=' ,E12.4,T50, 'N=',E12.4/ 

1 • NEWTON-RAPHSON METHOD : I T ER AT I ON #',13/ 

2 ' J' ,T15» 'XAB* ,T30, 'XAA* , T40,, s ' X AA (CALC) ',T55* 

3 *K( CALC ) ' ,T70, 'K( INPUT) ' ,T8 5 , • K-K (CALC ) • ) 

DO 700 J=1,JX 

IF ( XAAC! J ) .EO.XALF ) GO TO 600 
WRITE! I OUT, 7) J, XAB! J > , XAAC ( J ) , ALK ( J ) ,K0, R ( J > 

FORMAT! I5,E15.4,E30.4,3E15.4) 

GO .TO 700 

WRITE ( I OUT ,8)J,XAB( J ) , X AL F 

FORMAT! I 5,E15.4,20X, A4, ' ZERO NOT FOUND') 

XAAC! J)=-1.0 
CONTINUE 
WRITE! I OUT , 3 ) 

FORMAT! ' 1 ' ) 

CALL PLOT 1! NSC ALE, NHL ,NSBH ,NVL ,NSBV ) 

CALL PL0T2! GR ID, l.»0.»l.»0.) 

CALL PL0T3! •*' ,X AB , XAAC ♦ J X ) 

CALL PL0T4 ( 8, ; XAA' ) 

WRITE! I OUT , 2 ) 

FORMAT! /T25, ' XAB' ) 

WRITE! I OUT , 5 ) K0,M0,N0 

FORMAT!' PROGRAM REGSOL 2 ' / 1 K0= ' , E 1 3 • 5 , 10X » ' M0= ' » E 1 3. 5 * 1 OX , ' N0= ' , 
1 E 1 3 .5 ) 

RETURN 

END 



00016800 

00016900 

00017000 

00017100 

00017200 

00017300 

00017400 

00017500 

00017700 

00017800 

00017900 

00018000 

00018100 

00018200 

00018250 

00018300 

00018400 

00018500 

00018550 

00018600 

00018700 

00018800 

00018900 

00019000 

00019300 

00019400 

00019500 

00019600 

00019700 

00019800 

00019900 

00020000 

00020100 

00020200 

00020300 

00020310 

00020320 

00020330 

00020340 

00020400 

00020405 

00020410 

00020415 

00020420 

00020500 

00020600 

00020700 

00020840 

00020900 

00021000 

00021010 

00021020 

00020750 

00020760 

00020770 

00020780 

00020790 

00020800 

00020810 

00020820 

00020830 

00021300 

00021400 

00021450 

00021500 

00021600 

00021700 

00021800 

00021900 

00022000 

00022100 

00022200 

00022300 
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preceding page blank not eilkkd 

VI. PROGRAM MATRIX 
A. Purpose 


This is a general program to solve an equation of the type: 

y = aj X x + a 2 X 2 + . . . + a„ X n (9) 

(through the method of least squares) and, therefore, can be used for solutions 
of various problems. One example is in the solution of the following equation: 

taK sa = lnK D +4f<xJ - x“) +4f< 6 '*S *“ - ■> 

+ ^( X 5 " X i) + i7 (6 X A X B - D 

[5.5] 

which is an equation representing the distribution of a component between two 
asymmetric binary solutions. This program may also be used in place of 
REGSOL1. 


R. Numerical Method 


To solve (9) by method of least squares for n coefficients requires m > 
n observations of the form ( x ti , x 2i , . . . , x ni , y^ where 

a i x li + . . . + a n x ni = y i fore a chi = 1, 2, . . . , m. 

If A is an N x N matrix such that 

m 

A (i, j) = 2 x ik x ik , i = 1,2, ...,n, 

K - 1 y = 1, 2 , . . . ,n 

m 

and B(i) = 2 v k X ik i = 1 , . . . , n 

K = 1 

then Z = A -1 B, the solution of the matrix equation AZ = B, is the required least 
squares solution of (1) with Oj = Z{ , i = 1, . . . , n. This particular program 
allows n < 10, m < 50. 
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c. 


Notation used in Program MATRIX 


X 

X 

Y 

y 

NCOEF 

N 

NX 

M 

A 

A 

B 

B - before inversion of A 

B 

Z - following inversion of A 

Y0 

y - as calculated from (1) 

RY 

Y - Y0 


RSUM 



D. Input to and Output from Program MATRIX 

Card 1 : columns 1-5, 6-10 NCOEF, NX, respectively. 

Card 2 through NX + 1: columns 1-7, 8-14, . . . 71-77 X ki , k = 
1, . . . , NCOEF, Yj , respectively (fixed point format) . 

This set of cards constitutes a case. Cases may be stacked one behind the 
other. 

Output from each case is as follows: 

(1) I, (X(K,I), K = 1, NCOEF), Y(I) I = 1, 2, . . . , N 

(2) A -matrix, B-vector 

(3) Determinant of A, A -1 , Z 

(4) RSUM 

(5) (Y (I) , Y0(I), RY(I), I = 1, NX) 

Figure 9 shows a sample set of input while Figure 10 shows a sample set of 
output. 
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3 9 



1 . 

0.318 

0.958 

-3.183 

1. 

-.384 

.860 

-3.396 

1 . 

-.630 

.812 

-3.748 

1. 

-.728 

. 728 

-3.697 

1 . 

-.804 

.484 

-3.276 

1. 

-.798 

.318 

-2.845 

1. 

-.814 

-.066 

-2.145 

1. 

.050 

.934 

-3.278 

1 . 

-.078 

.942 

-3.667 


Figure 9. Sample Input to Program MATRIX 


X( 1 , * ) 

X< 2.*) 

<( 

3.* ) V 

l i.tff 

r . 31 8 


*955 - 3« 1 8 3 

2 l.cor- 

-O. 384 


• 86>t -3*3°e 

3 1 . * >K> 

-■'.63. 

ty 

• 812 -3*74 6 

a l . f 

- '. 72 8 ■ 


*728 — 3« b 9 7 

5 i • •: o 

-f . PI 4 


• 4 84 - 3» 2 7 e 

6 i • r <? '■ 

- ^ . 79 8 


.318 -2*945 

7 1. 

- s 31 4 

-O 

• ■!- 66 — 2m 14 5 

3 1.6 or- 

-.75'. 


• 9 34 - 3« 2 7 t 

9 i . : - 

-‘ .(78 


• <342 — 3© 66 7 

A( * , 1 ) 

A (* < 2 > 

4( * • 

3 > B 

0.950 f 1 — , 

o 36 70 1 

•5970 -1 - • 292D 2 

-0.337D CT 1 

.3130 .1 - 

• » 1 6BD 1 ^ • 12 ID i\ 2 

0.53 70 -1 — V 

a 1 6 80 " l 

>4550 1 - ~ • 2\ 60 2 

A-IN3F9SE NOW 

IN A# B CONTAINS THE! COEFFICIENTS 

D£TE3M= .6' 

U 90 C 1 



A( * , 1 ) 

A ( * * 2 ) 

A ( *, 

3 > rf 

C. 210 0 f- 1 ». 

• i5io i - 

.2 "3D 1 -'.'• 1540 el 

B • 1 5 1 0 1 

• i 4 8D •; i - 

.1320 .1 39 4D v* 

— 0, 2030 VI - 

• 13 2D • 1 

! .22 ' 

D '1 - . 1990 -..l 

SORT { 3ES ICUALS )= i. 296190 J 

. 

V 

YO 



-0.31 837 0 VI 

-0, 316710 

Cl -tf 

.158690- _*1 

-C.339 5C o r 1 

-1 . 35996D 

e i 1, 

.273560 7* 

-0.3743" C C,1 

372390 

» l -•> 

.241 430-31 

-0.3597' 0 Cl 

. 3644i 0 

c 1 -a 

.529 960-31 . 

-0.3276-0 71 

- .322540 

■*1 -t. 

,Si. 5 770— 31 

-0.23455C 1 

--. 2B89CD 

... 1 ( 

.440410- 11 

-0.2145(0 -1 

-t. 213760 

21 -7. 

.737 1 70— ? *2 

-0.3273CD < 1 

-r . 3 3597. 0 

7 1 C- 

. 819 € 40—11 

-0.356770 Cl 

. 348940 

<1 -C. 

. 1 77 6-ID Jfl 



AS FF 1 Y = B1 * XI +. .. +B-NC0EF4X-NC0EF 


Figure 10. Sample Output from Program MATRIX 



E. Listing of Program MATRIX 


c 

c 

c 

c 

c 


c 

100 

1 

8 

c 


,7 

2 

200 


,300 


500 

'600 

;s 

i 

1 800 
3 


6 


j 900 


1000 


11200 


A 


2000 


C 


PROGRAM MATRIX 
P.A.COMELLA 
CODE 641.1 

GODDARD SPACE FLIGHT CENTER 
GREENBELT .MARYLAND 20771 
IMPLICIT REAL*8 ( A-H.O-Z ) 

REAL*8 A ( 10,10) ,B( 10) 

1 »X(10,50)»Y(50),Y0(50)»RY(50) 

EQUIVALENCE (NCOEF.N) 

INTEGERS GAMMA! 10 ), ALPHA ( 10 ), BE TA , DELTA 

OATA GAMMA, ALPHA, BETA, DELTA/10*' X ( ' , 10* ' A ( * , ' , ' B • * • Y'/ 

INPUT 

READ! 5,1, END=2000) NCOEF.NX 
FORMAT (2 15) 

WR ITE! 6,8) ( GAMMA! I ) , I , 1=1 ,NCOEF ) .DELTA 
FORMAT! '1' , IX, A4, 12, ',*)', 10! 2X, A4, 12, ',*)') ) 

DO 200 1=1, NX 
INPUT 

READ (5,2) ( X! K, I ) ,K = 1 ,NCOEF) ,Y( I ) 

WRITE(6,7)I,!X(K,I) , K= 1 ,NCOEF ) , Y ( I ) 

FORMAT! I3,F8.3,10F11.3) 

FORMAT! 11F7.3) 

CONTINUE 
DO 300 K= 1 , N 
B ( K ) =0.000 
DO 300 KP= 1 ,N 
A! K,KP)=O.ODO 
DO 600 K= 1 , N 
DO 600 1 = 1, NX 
DO 500 K P= 1 , N 

A(KP,K)=A(KP,K)+X(KP,I )*X(K,I ) 

B(K)=B(K)+Y( I )*X(K, I ) 

WR I TE ! 6 , 5 ) (ALPHA! I), I , I =1 ,NCOEF ) .BETA 
FORMAT! 11 ( 3X,A4, 12, • )')) 

DO 800 K = 1 , N 

WR I TE ( 6, 3 ) (A(KP,K),KP=1,N),B(K) 

CONTINUE 
FORMAT! 11E11.3) 

CALL MAT1NV! A, N,B,1, DETERM) 

WR I TE ( 6, 6 ) DETERM 

FORMAT! 'OA-INVERSE NOW IN A.B CONTAINS THE COEFFICIENTS AS FF:', 

1 ' Y = Bl*Xl + ... + B-NCOEF*X-NCOEF'/' DE TE RM= ' , E 1 3 . 5 ) 

WR I TE ( 6 , 5 ) (ALPHA! I), I, I = 1 ,NCOEF ) .BETA 
DO 900 K = 1 , N 

WR I TE ( 6 , 3 ) (A(KP,K),KP=1,N),B(K) 

CONTINUE 

RSUM=0 .000 

DO 1200 1=1, NX 

SUM=O.ODO 

DO 1000 K=1,N 

SUM = SUM + B ( K ) * X ( K , I ) 

Y0( I )=SUM 

R Y ( I ) =Y< I ) — SUM 

R$UM = RSUM + RY< I )SRY< I ) 

CONT INUE 

RSUM=DSORT< RSUM ) 

WRITE (6, 4) RSUM,( Y( I ) , YO ( I ) , R Y ( I ) , I = 1 , NX ) 

FORMAT! '0SQRT(RESIDUALS) = ',E13. 5/' Y YO ' / 

1 (3E13.5)) 

GO TO 100 
RETURN 
END 

SUBROUTINE MAT INV! A ,N,B,M, DETERM) 

MAT I NV IS A VERSION OF THE SHARE SUBROUTINE OF THE SAME NAME 


00000300 

00000100 

00000300 

00000325, 

00000350 

00000400 

00000500 

00000525 

00000550 

00000600 

00000700 

00000725 

00000750 

00000800 

00000900 

00001000 

00001100 

00001200 

00001300 

00001400 

00001500 

00001600 

00001700 

00001800 

00001825 

00001850 

00001900 

00002000 

00002100 

00002200 

00002500 

00002525 

00002550 

O0OO2575 

00002590 

00002600 

00002700 

00002750 

00002800 

00002900 

00003000 

000031001 

00003200| 

00003300' 

00003400; 

00003500 

00003600; 

00003700 

000038001 

00003900; 

00004000; 

00004100 

00004200 

00004300 

00009700 
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IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 A( N,N> ,B( N,M) ,PI VOT( 10 ) 

I NT EGER*4 I P I VOT ( 1 0 ) , I NDEX < 10 , 2 ) 

EQUIVALENCE ( IROW,JROW) ,( I COL UM , JCOL UM ) , ( AM AX , T , S W AP ) 
DETERM= 1 .DO 
00 20 J=1,N 
20 I P I VOT ( J ) =0 

00 550 1=1, N 
AM AX=0 . 00 

DO 105 J=1,N 

IF( IPIVOTI JI.EO.l) GO TO 105 
DO 100 K= 1 , N 

I F ( IPIVOT(K)-l) 80,100,740 
80 IF(DA8S(AMAX).GE.0ABS( A( J,K) ) ) GO TO 100 
I ROW= J 
I COLUM=K 
AM AX=A( J , K ) 

100 CONTINUE 
105 CONTINUE 

1 PIVOT! I COLUM ) = I P I VOT ( ICOLUMl+l 
IF! IROW.EQ. ICOLUM) GO TO 260 
DETERM=-DETERM 

DO 200 L = 1,N 
SWAP=A ( I ROW , L ) 

A! I ROW * L ) = A ( ICOLUM, L ) 

200 A! ICOLUM, L)=SWAP 

IF(M.LE.O) GO TO 260 
DO 250 L=1,M 
SWAP = B ( IROW ,L ) 

B ( IROW , L ) = SW AP 
250 B< ICOLUM, L)=SV!AP 
260 INDEX! I, 1 ) = IROW 

INDEX! I , 2 ) = I COL UM 

PI VOT! I) =A( ICOLUM, ICOLUM) 

DETERM = OETERM*P IVOT! I ) 

A! ICOLUM, IC0LUM)=1. DO 
DO 350 L=1,N 

350 A! ICOLUM, L)=A( I COLUM , L ) / P I VOT ( I ) 

IF(M.LE.O) GO TO 380 
DO 370 L = 1 , M 

B ( ICOLUM, L)=B( ICOLUM,L (/PIVOT! I ) 

370 CONTINUE 
380 DO 550 L 1 = 1 ,N 

IF! LI. EO. ICOLUM) GO TO 550 
T = A(L1, ICOLUM) 

A! LI, ICOLUM ) =0 . DO 
DO 450 L = 1 , N 

450 A(L1,L)=A(L1,L)-A( ICOLUM,L)*T 

IF(M.LE.O) GO TO 550 
DO 500 L = 1 , M 

500 B ( L 1 , L ) =B ( L 1 , L ) -B( ICOLUM,L)*T 
550 CONTINUE 

DO 710 1=1, N 
L=N+ 1- I 

IF! INDEX! L,1 ) .EO. INDEX! L ,2) ) GO TO 710 
JROW= I NDE X ( L , 1 ) 

JCOLUM= INDEX! L ,2 ) 

DO 705 K=1,N 
SWAP=A ( K , JROW ) 

A(K,JROW)=A(K,JCOLUM) 

A(K,JCOLUM)=SWAP 
705 CONTINUE 
710 CONTINUE 
740 RETURN 
END 


00009800 

00009900 

00010000 

00010100 

00010200 

00010300 

00010400 

00010500 

00010600 

00010700 

00010800 

00010900 

00011000 

00011100 

00011200 

00011300 

00011400 

00011500 

00011600 

00011700 

00011800 

00011900 

00012000 

00012100 

00012200 

00012300 

00012400 

00012500 

00012600 

00012700 

00012800 

00012900 

00013000 

00013100 

00013200 

00013300 

00013400 

00013500 

00013600 

00013700 

00013800 

00013900 

00014000 

00014100 

00014200 

00014300 

00014400 

00014500 

00014600 

00014700 

00014800 

00014900 

00015000 

00015100 

00015200 

00015300 

00015400 

00015500 

00015600 

00015700 

00015800 

00015900 

00016000 

00016100 

00016200 
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VII. PROGRAM QUASI 
A. Purpose 


If we have data on a complete distribution isotherm, this program may 
be used to find 2W/ZRT for each of the two crystalline solutions assuming that 
they are regular solutions with the quasi-chemical approximation. 


B. Numerical Method 


The basic equation is 


K = K r 


(1 + * 1A tf- D) 

Z 1^2 

[ 1 + 0 2B (P ~ 1>V 

2 ! 

i ^2 A » + ') ] 


( *)) 

(1 +* 2A (P- D] 

i¥- ) 

1 Mi.U'- i>l 

l *, A » + i) ! 

( 1 

l * 2B W + 1) ) 


z 1 qJ 


z 2l2 


( 10 ) 

where the symbols 0 1A , 0 2A , 0 1B , 0 2B , 0, 0', q 1; q 2 , q' x , q 2 , Z x and Z 2 cor- 
respond to 0“, 0g, 0^, 0g, 0^, q“, q“, q^, q|, Z a and respectively in 

Saxena's equation [5.6]. Letting 

Zq x 
Aa 


'2A 


Ab 



; 2B 


1 + 0 1B (fi'~ 1) 
*2B W + 1) 


-i zq ’ 2 
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Equation (10) may be written as 


K = 

which is the notation used in the 


^1A X 2B f l A < *2B 
^2 A ^1B ^2 A ^1 B 

program. Now 


H> + 4X 1A X 2A ( e ^-')} B 

0'={> - 4X 1B X 2B (e^ F - >)}* 


The problem is to find 


2c o Zoo 

ZRT ’ ZRT ' 


HD 


( 12 ) 


This is done by the method of non-linear least squares as outlined in the follow- 
ing paragraphs : Write (11) 


/( x i A , x iB> k ',y,y') = k' 


v f f 

^2B 1 1 A ^ l 2 B 

^1B *2 A f lB 


( 12 ) 


where 


Set 


Let 


_ J_ _ 2oj , _ 2o>' 
K ’ y " ZRT ’ y “ ZRT 


fobs ^ 2 a/X 


1 A > 


V = /( X 1A , X 1B , k'.y.y') - f obs 


(13) 


k’ - kJ'Q-j + Ak’ 

y = y( 0 ) + A y 
y' = y'(o) + A y' 

where k (o)> y(o)> y’(o>are initial estimates of k', y, y'. 
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Then linearize (13) by doing a Taylor Series expansion about 

( k (o)>y(o).y’(o)) = (°) 

so obtaining 


v + fobs - /( x iA»XiB» k '(o)>y(o)»y'(o)) 




(o) 


+ 5y — , 

(o) fi y 


( 14 ) 


(o) 


and now use method of least squares to solve this equation, iterating until 

|5k'|, |5y|, |5y'| 


are less than some prescribed e . 


C . Notation used in Program QUASI 


Z 1 5 

Z1 



z 2 : 

Z2 



q, : 

Q1 

q 2 s 

Q2 

q', : Q1P 

q ’ 2 : 

Q2P 

k (o) : 

KO 

y(o) : YO 

of observations of the pairs : 

(XiA* X 1B ): NX 

X 1A 

X1A 

X 2A 

X2A 

01A 

PHI1A 

02 A 

PHI 2 A 

01B 

PHI IB 

02 B 

PHI2B 

P 

BETA 

P 

BETAP 

dl 

DBDY 

d<3 

DBDYP 

dy 


dy' 


flA 

FlA 

flA 

F2A 

f IB 

FIB 

^2B 

F2B 

d/iA 


d/2 a 


dy 

DF1ADY 

dy 

DF2ADY 

flA 





FA12 



f2A 




d/j b 


d/ 2B 



DF1BDY 


DF2BDY 

dy 


dy' 



y(o) 


YPO 
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FBI 2 


hz 

/lB 

JL 

3k' 

3y 
3y' 

/ (X 1A , X 1B , k', y, y’) 
fobs • ^ 


DFDK 


e- : DFDY 


DFDYP 


F 


Determinant of Least Square Matrix : DET 


Least square coefficients 


Solution: 


6k' 

: Al 

k' 

: KPO 

k 

: KP 

y 

: YOO 

y' 

: YPOO 


5y : A2 


5y' : A3 


D. Input to and Output from Program QUASI 


Input: 

Card 1 : columns 1-5, 6-10, . . . , 41-45 Zl, Z2, Ql, Q2, QlP, Q2P, 
KO, YO, YPO, respectively (fixed point format). 

Card 2: : columns 1-5 NX 


Cards 3 & ff : columns 1-5, ... , 76-80 (X 1A d>, X^W), 1 = 1, 2, 
. . . , NX fixed point format, 8 pairs (X 1A , X 1B ) per card 
until NX pairs entered with a maximum of 50 pairs. 

Output from Program QUASI as follows : 

(1) The input 

(2) for each iteration: KP, KPO, YOO, YPOO, DET, Al, A2, A3 
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(3) for the final iteration 

X 1A (input), X 1A (calculated), / 1A , / 2A , p (y) 

X 1B (input), X 1B (calculated) » /ib ./2b » P ' (y) 

RE SID : Y-F 

Figure 11 shows a sample set of input while Figure 12 shows a sampling of 
output. 


2.0 2.0 1.0 1.0 1.0 1.0 0.15 1.5 1.0 

6 

0 . 0210 . 3410 . 0700 . 6920 . 1110 . 7990 . 1500 . 8500 . 2620 . 8970 . 3370.902 

Figure 11. Sample Input Data for Program QUASI 
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CONSTANTS AND INITIAL VALUta 


Z 1 = 

0.20000E 

01 

Z2= 

0.20000E 

0 1 

Ql — 

0. 1 OOOOE 

01 

02= 0.1 OOOOE 01 

01P= 

0. 1 OOOOE 01 

K0 = 

0. 15000E 

00 

Y0 = 

0. 1S000E 

01 

YPO= 

0. 1 OOOOE 

01 




I 



XI A 



XI B 






1 



0.021 



0.341 






2 



0.070 



0.692 






3 



0.111 



0.799 






4 



0.150 



0.850 






5 



0.262 



0.897 






6 



0.337 



0.902 







iteration: 

0 

K0 = 

0. 

1 50 0 OE 

00 

1/K0= 

0.66667E 

01 

Y= 

0. 15000E 

01 

YP= 

0. 10000E 

01 



OE T = 













i terat ion: 

1 

K0 = 

0. 

1084 IE 

00 

l/K0 = 

0. 9224 1 E 

01 

Y= 

0. 1 1791E 

01 

YP= 

0.61774E 

00 



DET= 

0. 

1 1663E 

04 

A 1 = 

0.25574E 

01 

A 2= 

-0.32086E 

00 

A 3= 

-0.38226E 

00 

I teration: 

2 

KC = 

0. 

1 07 6 1 E 

00 

l/KO = 

0.92928E 

01 

Y= 

0. 1 2412E 

0 1 

YP = 

0 .70080E 

00 



DET = 

0. 

19199E 

04 

A 1 = 

O.68707E- 

01 

A 2= 

0.6 2034 E 

-0 1 

A 3= 

0.830S6E- 

01 

i teration: 

3 

K0 = 

c. 

10806E 

00 

1/K0 = 

0.92537E 

01 

Y= 

0. 1 2450 E 

01 

YP = 

0.701 17E 

00 



DET = 

0 • 

23756E 

CA 

A l = 

— 0 • 391 52E— 0 1 

A 2= 

0.37947E 

-0 2 

A 3= 

0.37367E- 

03 

i teration: 

4 

K0 = 

0. 

1080.3E 

CO 

1/K0 = 

0.92563E 

01 

Y= 

0. 12446E 

01 

YP = 

0.70100E 

00 



OET = 

0. 

23633E 

04 

A 1 = 

0.2S766E- 

02 

A 2 — 

-0.35054E-03 

A3= 

— 0 . 17069E— 

03 

i teration: 

5 

K0 = 

0. 

10804E 

00 

1/K0= 

0 . 9256 1 E 

01 

Y = 

0. 12446E 

01 

YP= 

0 • 70 1 0 1 E 

00 



DET = 

0. 

236 15E 

04 

A 1 = 

-0.12327E- 

03 

A 2= 

0. 17373E-04 

A 3= 

0.77288E- 

05 

i teration: 

fc 

K0 = 

0. 

10804E 

00 

l /K0= 

0.92562E 

01 

Y = 

0. 12446E 

01 

YP = 

0.70 10 IE 

00 



GET = 

0 . 

23569E 

04 

A 1 = 

0. 81 759E- 

05 

A 2= 

0 . 65433E— 06 

A 3= 

-0.25030E- 

05 


G2P= 0.1 OOOOE 01 


FINAL RESULTS 
X1A-INPUT 
0 .21000E-01 
0.70000E-01 
O.llIOOE 00 
0*1 5 OOOE 00 
0.26200E 00 
0.33700E 00 


XI A-CALC 
0.2 101 IE- 0 1 
0.665 1 IE-01 
0.1 139 3E 00 
0.1 5956F. 00 
0.26104E 00 
0.30166E 00 


F 1A 

0.31550E 01 
0.26420E 01 
0.23503E 01 
0.21 3976 01 
0.17J52E 01 
0.15578E 01 


F2A 

0. 10010E 01 
0. 10093E 01 
0.1 02 1 OE 01 
0. 10355E 01 
0.10927E 01 
0 • 1 1 44 1 E 01 


0ETA( Y) 
0.10969E 01 
0.12820E 01 
0.14056E 01 
0.15035E 01 
0.17064E 01 
0.17914E 01 


X IB- INPUT 
0.34100E 00 
0.69200E 00 
0 . 7990 OE 00 
0 . 85000 E 00 
0.89700E 00 
0.90 200 E 00 


FIB 

0.13107E 01 
0.10689E 01 
0.10314E 01 
0.1018 3E. 01 
0.10091 E 01 
0.10083E 01 


F2B 

0.10832E 01 
0.13476E 01 
0.14966E 01 
0.15892E 01 
0.16925E 01 
0.17047E 01 


BETAPI YP» 

O. 1383 IE 01 
0.1 366 OE 01 
0.128S5E 01 
0. 12321E 01 l 
0.11728E 01 
0.116S8E 01 


RESID 

-0.10975E-04 
0.1 4887E— 02 
—0. 29345 E— 02 
“0. 95568 E— 02 
0.96017E-03 
0. 35336E— 0 l 


Figure 12. Quasi-Chemical Approximation: Equation (5.6) 
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E. Listing of Program QUASI 


c 

C PROGRAM QUASI 

C QUASI CHEMICAL APPROXIMATION 

C P.A. COMELLA 

C CODE 641.1 

GODDARD SPACE FLIGHT CENTER 
GREENBELT .MARYLAND 20771 


50 

1 

3 


5 


200 

900 


REAL *4 KO.KP.KPO.XIA! 50 > ,X1B< 50) ,PHI 1B( 50) .PHI 2B< 50) ,X2B< 50) . 00000100 

1 X21BI 50) ,X1AC< 50 ) ,FA1 ( 50 ). FA2 (50) ,FB1 ( 50) ,FB2< 50) ,BET( 50) , 00000200 

2 BETPt 50 ) ,FUN( 50 ) ,DI FF( 50) .RES1DI 50) 00000300 

3 . X B 1 2 ( 50) 00000325 

LOGICAL* 1 TEST 00000350 

INTEGER *4 IN/5/. OUT/6/ 00000375 

COMMON SUM 1 , SUM 2 , SUM3 , SUM4 , SUM 5 , SUM6 , SUM 7 , SUMS . SUMO 00000400 

READ! IN, 1 , END = 2000 > Z 1 , Z2 ,01 ,02, 01P ,Q2P, KO, YO, YPO, NX , < X 1A (I) , X IB < I 100000500 

1 ,1 = 1, NX) 00000600 

FORMAT! 9F5.0/I5/I 16F5.0) ) 00000700 

WRITE ( OUT , 3 ) Z1 , Z2 , 0 1 ,02 ,01 P ,02P , KO , YO, YPO , ( I , XI A! I ) , 00000710 

1 X 1 B ( I ) , I = 1 , i'J X ) 00000720 

FORMAT ( 1 1 1 , 10X , ' OUAS I -CHEM I CAL APPROXIMATION: EOUAT I ON (5. 6)'/ 00000730 

1 'OCONSTANTS AND INITIAL VALUES'/ 00000740 

2 1 71= 1 ,E13.5,5X, 'Z2= ' ,E13.5, 5X, *01 = ' ,E13.5, 5X, '02= 1 , 00000750 

3 E13.5.5X, ' 01P=' , E13.5, 5X, '02P=' ,E13.5/ ' K0= • ,E 13. 5, 5X, ■ Y0= ' , 00000760 

3 E13.5.4X, 1 Y P 0 = ' ,E13.5/'0 I',15X, 00000765 

4 ' X1A 1 , 17X, ' X1B' /( I 5, 10X.F10.3, 10X, F10.3 ) ) 00000770 

712=0.5*71 00000800 

Z 22=0. 5*72 00000850 

Z01 =Z 12*01 00000900 

Z 02 = 7 2 2*02 00001000 

ZOP 1 = Z 1 2 *0 1 P 00001100 

Z0P2 = 7 22*02 P 00001200 

Z011=Z01-1.0 00001300 

Z 0 2 1 = 7 0 2 - 1 . 0 00001400 

Z Q P 1 1 = Z Q P 1 - 1 . 0 00001 500 


Z 0 P 2 1 = 7 Q P 2 - 1 • 0 
K P = 1 . 0 / K 0 
KP() = KP 
Y 0 0 = Y 0 
YPOO= YPO 
R S U M = 0.0 
TEST=. FALSE. 

I NDEX = 0 
WRITE (OUT, 5) 

FORMAT! '-' ) 

WRITE! (JUT , 4 ) INDEX, KO.KPO.Y 00, YPOO 
I ND EX = 1 
DU 200 1=1, MX 
X? B ( I ) = 1 .O-Xlh ( I ) 

PHI IB! I ) = X 1 B ( I ) * 0 1 P / ( X 1 B ( I ) * 0 1 P + X 2 B ( I )*U2P ) 
PHI2B! I >=1. O-PHI ltf! I ) 

X21B! I )=X2B( I ) / X 1 K ( I ) 

X 1 AC ( I ) = X 1 A ( I ) 

X 8 1 2 ( I ) = 2 . 0*X 1 B ( ! ) * X 2 B ( I ) 

CONTINUE 
SUM 1 = 0.0 
SUM 2 = 0 . 0 
SUM 3 = 0.0 
SUM4 = 0 .0 
SUM 5 = 0.0 
SUM 6= 0.0 
SUM 7 =0.0 
SUM 8 = 0.0 
SUM9=0. 0 
fc Y = F X P ( YOU ) 


00001600 
00001700 
00001800 
00001900 
00002000 
00002025 
00002050 
00002075 
00002080 
00002085 
00002090 
00002095 
000.0 2100 
00002200 
00002300 
00002400 
00002500 
00002600 
00002650 
00002700 
00003000 
00003100 
00003200 
00003300 
00003400 
00003500 
00003600 
00003700 
00003800 
00003810 
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E y P = E X P ( Y POO ) 

EY1=EY— 1.0 
EYP1=EYP-1.0 
DO 1000 1=1, NX 
XA1 = X1A( I ) 

XA2 = 1 .0-XA1 

XA12=2.0*XA1*XA2 

BET A=SQRT ( 1 .0+2 .0*X A 1 2*E Y 1 ) 

BP 1 = 6E T A + 1.0 
BM1=BETA-1 .0 

BETAP = SORT( 1. 0 + 2.0 YXB12I I )*EYP1 ) 

PHI lA = XAl#01/( XA1*01+XA2*U2 ) 

PHI2A=1. O-PHI 1A 
BPP1=BETAP+1 .0 
BPM1 = BET AP-1 . 0 
BPP10=BPP1*BPP 1 
BP10=BP1*BP1 
DBDY=XA12*E Y/BETA 
DH0YP = XB 1 2 ( I )*EYP/BETAP 
FlA=1.0+PHI2A*BHl/( PHI 1A-BP1 ) 

F2A=1.0+PHI lA*BMl/< PHI2A*BP1 ) 

F1B = 1 .0+PH12BI I >*BPMl/( PHI 1B( I )*BPP1 ) 

F2B=1.0+PHI 1 B ( I >*BPMl/( PHI2BI I )*BPP1 ) 

DF 1 AD Y = Z01*F 1 A**ZO 1 1 *2 .0*PH 1 2A / ( PH I 1A*BP10 )*DBDY 

DF2A0Y=Z02*F2A**Z021*2.0*PHI 1A/ ( PH 1 2 A*BP lu ) *DBD Y 

F1A=F1A**ZQ1 

F2A=F2A**Z01 

FA12=F1A/F2A 

DF1BDY=Z0P1*F1B**Z0P11*PHI2B( I I/IPHI 1 H < I ) *BPP 1 0 ) *DriO YP 
1*2 .0 

DF2BDY=ZOP2*F2B**Z(OP21*PHI 1 B ( I >/<PHI2H( I ) *BPP 1 0 ) *OnU YP 
1*2.0 

F1B=F1B**ZQP1 
F2B=F2B**ZQP2 
FB12=F2B/F1B 
•DFDK = X2 IB ( I)*FA12*FB12 

DFDY = KP0*X21B( I ) *FB 1 2* ( F2A*DF 1 ADY-F 1 A*0F2 ADY ) / ! F2 A = F?A ) 
DFDYP=KP0*X21B( I ) * F A 1 2* ( F 1 B*D F 2BD Y- F 2B*DF 1 BD Y ) / ( F 1 h * F 1 B ) 
F=KPO*DFDK 
Y = XA2/X1A ( I ) 

YF=Y-F 

X 1 AC ( I ) =1 .0/ ( 1.0+F) 

IF (TEST .EQ. .FALSE.) GO TO 950 

F A 1 ( I ) = F 1 A 

FA2< I ) = F 2 A 

FBI ( I ) = F1B 

FB2 ( I ) =F 1 2B 

FB2 ( I ) =F2B 

BET ( I ) =BE T A 

BET P ( I ) =BETAP 

FUN ( I )=F 

DIFFI I )=YF 

RES I 0 ( I ) = X 1 A ( I ) -X 1 AC ( I ) 

RSUM = RSUM+RESID( I ) *RES I D ( I ) 

C GO TO 1000 

950 CONTINUE 

SUM 1=SUM1+DFDK*DF0K 

SUM2=SUM2+0FDK*DFDY 

SUM3=SUM3+DFDK*DFDYP 

SUM4=SUM4+DFDK*YF 

SUM5=SUM5+DFDY*DFDY 

SUM6= SUM6+0FD Y*0 FD Y P 

SUM7=SUM7+DFDY*YF 

SUM8=SUM8+DFDYP*DFDYP 

SUM9=SUM9+0FDYP* YF 


00003820 
00003830 
00003890 
00013900 
00014000 
00014100 
00014200 
000.14400 
00014500 
00014600 
00014700 
00014725 
00014750 
00014800 
00014900 
00015000 
0001 5100 
00015300 
00015400 
00015500 
00015600 
00015700 
00015800 
00016100 
00016200 
00016300 
00016400 
00016500 
00016600 
00016650 
00016700 
00016715 
00016725 
00016750 
00016800 
00016900 
00017000 
00017100 
00017200 
00017225 
00017250 
00017300 
00017400 
00017500 
00017600 
00017700 
00017800 
00017900 
00018000 
00018100 
00018200 
00018300 
00018400 
00018500 
00018600 
00018700 
00019000 
00019100 
00019200 
00019300 
00019400 
00019500 
00019600 
00019700 
00019800 
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1000 CUNT I NUE 

IF (TEST ,E0. .TKUE.) GO TO 1100 

DET = DETEKM( 1 ,5, 8,6,6 )+DETERM< 2.6, 3,2,8 >+DETERM< 3,2,6, 5,3) 
Al=(DETERM(4,5,8,6,fe)+DETERM(7,6,3,2,8)+UETERM(9,2,6,5,3))/DET 
A2=( DETERMI 1 , 7 , 8 ,6 ,9 ) +DETERM ! 2 .9*3* A, 8 1+DETERMl 3,4,6, 7,3 ) ) /GET 
A3 = ( DETERM 1, 5,9, 6,7) + DE TERM! 2, 4, 6, 2, .9 )+DETERM( 3,2,7,5,41 ) /UET 
TEST=( AHS( Al/KRO ) .LT.l .E-05) .AND. ( AttSf A2/Y00 ) .IT. 1 .E-Ob ) .AND. 

1 ( ABS( A3/YP00 ) .IT. 1 .E-05 ) .UR. INDEX.GT. lo 

K PO = K P0 + A 1 
Y00= Y00+ A2 
Y P 00 = Y P 00 + A 3 
KP=1.0/KPO 

WHITE! OUT ,4 ) INDEX ,KP,KPO, 

1 YOO,YPOO,DET,A1,A2,A3 

I F( TEST .EO. .TRUE . ) GO TO 900 
I I'lDE X = I NDE X+ 1 

A FORMAT ( 1 ITERATION: ', 13, 5X, 'K0=' ,E13.5,5X, ' 1/K0= ' ,E13.5, 

1 5X,'Y=',E13.5,5X,'YP=',E13.5/18X,'0ET=',E13.5,5X, 

2 ' A 1 = • ,E13.5,4X, 'A2=' ,E13.5,5X, , A3 = ',E13.5) 

GO TO 900 

1100 KP=1.0/K PO 

WRITE ( OUT, 6) ( XI A ( I ) , XI AC( I ) ,FA1 ( I ) , FA? ( I ) , BET ( I ) , 

1 X 1 B ( I ) »FB1( I ) »FB2( I ) »BETP( I )»KESID( I ) ,1 = 1* NX) 

6 FORMAT! '-FINAL RESULTS'/' X 1 A- I NPUT ' , 4X , ' X 1 A-C ALC ' , HX , 

1 ' F1A ' , 8X, ' F2A ' ,6X , ' BETA ( Y ) ' ,3X , ■ X lb- INPUT ' , 7X, 

2 ' FIB' ,8X, ' F2B ' ,5X, 'BETAP(YP) ' ,2X, ' RE S I D ' / 

3 ( 10E13.5) ) 

GO TO 50 

2000 RETURN 
END 

REAL FUNCTION DETERM *4( I ,J,K,L,M) 

COMMON SUMI20) 

DET ERM = SUM ( I ) *( SUM( J ) *SUM ( K ) -SUM ( L )' * SUM ( M ) ) 

RETURN 

END 


00019900 

00019950 

00020000 

00020100 

00020200 

00020300 

00020400 

00020500 

00020700 

00020800 

00020900 

00021000 

00021015 

00021025 

00021040 

00021050 

00021075 

00021078 

00021081 

00021100 

00021200 

00021300 

00021400 

00021500 

00021600 

00021625 

00021650 

00021700 

00021800 

00021900 
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